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We overview a series of recent works devoted 
to variance reduction techniques for numerical 
stochastic homogenization. Numerical homogenization 
requires solving a set of problems at the micro 
scale, the so-called corrector problems. In a random 
environment, these problems are stochastic and 
therefore need to be repeatedly solved, for several 
configurations of the medium considered. An empirical 
average over all configurations is then performed 
using the Monte-Carlo approach, so as to approximate 
the effective coefficients necessary to determine the 
macroscopic behavior. Variance severely affects the 
accuracy and the cost of such computations. Variance 
reduction approaches, borrowed from other contexts 
of the engineering sciences, can be useful. Some of 
these variance reduction techniques are presented, 
studied and tested here. 
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We overview a series of recent works related to some random multiscale problems motivated by 
practical problems in Mechanics. For simplicity, we argue on a linear elliptic scalar equation in 
divergence form: 

I -div[A(J)Vn^]=/ in V, 

I = 0 on dV, 

although the scope of the techniques we describe go beyond this simple setting. The matrix 
coefficient A is assumed random stationary. The purpose is to compute the homogenized matrix 
coefficient A* present in the homogenized equation 


-div =f in V, 

u* = 0 on dT>, 


( 0 . 2 ) 


which captures the average behavior of the solution to (0.1). 

We begin by recalling in Section 1 the basics of homogenization theory, both in the 
deterministic (periodic) context and in the random context, which are useful for our exposition. 
Next, we successively present three different variance reduction techniques for the problem 
considered. We begin in Section 2 with the classical, general purpose technique of antithetic 
variables. The efficiency of that technique is substantial, but is also limited in particular because 
the technique does not exploit much the specifics of the problem considered. We present in 
Section 3 the technique of control variate, which requires a better knowledge of the problem 
at hand. A problem simpler to simulate and close to the original problem, in a sense that is 
made precise below, has to be considered and concurrently solved. The technique uses that 
knowledge to, effectively, get a much better reduction of the variance. In Section 4, we expose 
a slightly different approach, imported from solid state physics, namely that of special qiiasi- 
random structures. It consists in selecting (somewhat in the spirit of another well-known technique, 
stratified sampling) some configurations of the random environment that are more suitable than 
generic configurations to compute the empirical averages, so as to again minimize the variance. 
Our final Section 5 presents some further research directions. 


Before we proceed, we mention that we will assume throughout our text that the reader is 
reasonably familiar with the homogenization theory. We refer to the classical textbooks [5,12] 
for this topic. We also mention [1,13,14] for general presentations and overviews of the issues 
examined here, along with some related issues. 


1. Brief overview of classical homogenization settings 


(a) Periodic homogenization 

To begin with, we recall some well known, basic ingredients of elliptic homogenization theory 
in the periodic setting. We consider (0.1) in a regular, bounded domain V in R'^, and choose the 
matrix coefficient A = Aper to be symmetric and Z'^-periodic. Note that, throughout our text, we 
manipulate for simplicity symmetric matrices, but our discussion in Sections 2 through 4 carries 
over to non symmetric matrices up to slight modifications. 

The corrector problem associated to (0.1) in the periodic case reads, for p fixed in R'^, 


-div (Aper (p) (p -b Vwp)) = 0, 
Wp is Z'^-periodic. 


( 1 . 1 ) 


It has a unique solution up to the addition of a constant. Then, the homogenized coefficient in (0.2) 
reads 


= 


ef Aperiy) {cj + Vwej (p)) dy, 
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where Q = (0,1)'^ is the unit cube and (ei)i<i<d is the canonical basis of Equivalently, A* 
satisfies 


VpeR"*, A*p 


Aperiv) {p + Vwp(i/)) dy. 

Q 


The main result of periodic homogenization theory is that, as e goes to zero, the solution u® to (0.1) 
converges to u* solution to (0.2). The convergence holds in lF'{T>) and weakly in Hq(T>). The 
correctors may also be used to "correct" m* in order to identify the behavior of in the 
strong topology of Hq {T>) . 


Practically, at the price of only computing d periodic problems (1.1), the solution to 
problem (0.1) can therefore be efficiently approached for e small. 


(b) Stochastic homogenization 

Because this is well known and for the sake of brevity, we skip all technicalities related to the 
definition of the probabilistic setting (we refer e.g. to [1] for all details). We assume that A is 
stationary in the sense 

VfceZ'^, A{x + k,Lo) = A{x,t]^u)) almost everywhere in X, almost surely (1.2) 


(where r is an ergodic group action). We consider the boundary value problem (0.1) for 
A — A{-,ui). Standard results of stochastic homogenization [5,12] apply and allow to find the 
homogenized problem. These results generalize the periodic results recalled in Section (a). The 
solution M®(', uj) to (0.1) converges to the solution to (0.2) where the homogenized matrix is now 
defined as 


[A*]ij = E ^ A{y,-) [cj + Vwe^- {y, ■)) dyj , 


where, for any p G R"^, Wp is the solution (unique up to the addition of a random constant) to 


—div [,4 (i/,a;) (p + Vtyp(j/,tj))] = 0 a.s. onR'^, 
Vwp is stationary in the sense of (1.2), 


E 




(1.3) 


Note that is a random function, while its homogenized limit u* is deterministic since A* is 
deterministic. 

A striking difference between the stochastic setting and the periodic setting can be observed 
comparing (1.1) and (1.3). In the periodic case, the corrector problem is posed on a bounded 
domain (namely, the periodic cell Q), since the corrector Wp is periodic. In sharp contrast, the 
corrector problem (1.3) of the random case is posed on the whole space R"^, and cannot be reduced 
to a problem posed on a bounded domain. The fact that the random corrector problem is posed on 
the entire space has far reaching consequences for numerical practice. Truncations of problem (1.3) 
have to be considered. The actual homogenized coefficients are only captured in the asymptotic 
regime. 

More precisely, the deterministic matrix A* is usually approximated by the random matrix 
A%{lo) defined by 


Vp G ] 


A^(a;)p = 


IQatI 


A{-,uj) (p + Vwp 


(1.4) 


which is obtained by solving the corrector problem on a truncated domain, say the cube Qn = 

— div |^A(-,a;) ^p + Vro^(-, o;)^ j = 0, Wp {-jUj) is Qn- periodic. (1.5) 

Although A* itself is a deterministic quantity, its practical approximation A^ is random. It is 
only in the limit of infinitely large domains Qn that the deterministic value is attained. As shown 
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in [8], we have 


lim A%{ui) = A* almost surely. 

N —¥00 


As usual in the random context, the error A* — A^{uj) may be expanded as 

A* - A*^{co) = (a* - E [A^] ) + (e [A^] - A^(o.)), (1.6) 

that is the sum of a systematic error and of a statistical error (the first and second terms in the above 
right-hand side, respectively). 

A standard technique to compute an approximation of E [A^] is to consider M independent 
and identically distributed realizations of the field A, solve for each of them the corrector 
problem (1.5) (thereby obtaining from (1.4) i.i.d. realizations A^™(a;), 1 < m < M) and compute 
the Monte Carlo approximation 

M 

m—1 

In view of the Central Limit Theorem, we know that E asymptotically lies within the 

confidence interval 


rMC 

JM 


1.96- 




tNIC I 

,Im -1-1.96- 


y/M 


with a probability equal to 95 %. 

For simplicity, and because this is overwhelmingly the case in the numerical practice, we 
have considered in (1.5) periodic boundary conditions. These will be the conditions we adopt 
throughout our study. Other boundary conditions, or approximations, may be employed. The 
specific choice of approximation technique is motivated by considerations about the decrease of 
the systematic error in (1.6). Several recent mathematical studies by A. Gloria and R Otto [11] have 
clarified this issue. The variance reduction techniques we present in this article can be applied to 
all types of boundary conditions. 


2. Variance reduction using antithetic variables 

We present here a first attempt [6,7,9] to reduce the variance in stochastic homogenization. The 
technique used for variance reduction is that of antithetic variables. 

The variance reduction technique using antithetic variables consists in concurrently 
considering two sets of configurations for the random material instead of only one set. The two 
sets of configurations will be deduced one from the other. Indeed, fix M — 2A4. Suppose that we 
give ourselves A4 i.i.d. copies of A(x,lu). Construct next A4 i.i.d. antithetic 

random fields 

B™'{x,lo) = T (^A^{x,lo)) , l<m<M, 

from the (A'"(a;,a;))j^<^<^. The map T transforms the random field A™ into another, so-called 
antithetic, field S'". The transformation is performed in such a way that, for each m, has the 
same law as A™, namely the law of the matrix A. Somewhat vaguely stated, if A was obtained in 
a coin tossing game (using a fair coin), then B"* would be head each time A"* is tail and vice versa. 
Then, for each 1 < m < A4, we solve two corrector problems. One is associated to the original field 
A™, the other one is associated to the antithetic field B"*. Using its solution we define the 
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antithetic homogenized matrix the elements of which read, for any l<i,j<d, 




IQtvI 

And we finally set, for any 1 < m < M, 


T rjTn 

e,- B 


(•,^) (- 


Cj + Vw, 


N,m 


(■.^)) ■ 


Since A™ and are identically distributed, so are and B^™. Thus, A^"* is unbiased (that 
is, E = E (A^™)). In addition, it satisfies: 

A*i™ —> A* almost surely, 


N 


N—^-\-oo 


because A™ and are ergodic. The hope is that the new approximation A^"* has less variance 
than the original one A^"*. It is indeed the case under appropriate assumptions. 

The approach has been studied theoretically in [6,7,9], in the one-dimensional setting and 
in some specific higher dimensional cases. The approach is shown to qualitatively reduce the 
variance. A quantitative assessment of the reduction is however out of reach. Only numerical 
tests can provide some information in this direction. 


The tests we have performed in [6,9] concern various "input" random fields A{-,uj), some 
i.i.d., some correlated, with various correlation lengths. In these settings, we have investigated 
variance reduction on a typical diagonal [A]^(a;)]]^]^, or off-diagonal [A]^(a ;)]]^2 entry of the 
approximate homogenized matrix A^(a;), as well as on the eigenvalues of the matrix A^{uj), 
and the eigenvalues of the associated differential operator L = — div [A]^(a;)V-] (supplied with 
homogeneous Dirichlet boundary conditions on dV). 

Let us give one such example. Consider, in dimension two, the matrix A{x, cu) defined by 


A(x,uj) 


I Q 

feez2 



( 2 . 1 ) 


where Q — (0,1)^, (afe)fcgz 2 is an i.i.d sequence of random Bernoulli variables such that P(a/^ = 
a) = P(afc = /3) = 1/2, with a = 3 and /3 = 20. An example of the realization of each matrix field 
A{x, uj) and B(x, uj) is given in Figure 1 (in black, the value a and in pink, the value /3). 



Figure 1. An example of realization of A(x, uj) together with its antithetic field B{x, uj) (reproduced from [9]). 
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We then compare two computations with identical cost. For this purpose, we first use a 
classical Monte Carlo method with 2M draws (with here 2M = 100). Second, we apply the 
antithetic variable technique using only M draws. Since we solve two corrector problems for 
each of the draws (one for A™ and one for S'"), the numerical cost is equal to the cost of the 
classical computation. The results are shown in Figure 2, where we can see that the (numerically 
estimated) variance is reduced. 



Size of 


Figure 2. Estimation of A*^ (with confidence interval) with respect to |Q]v| (in red, the classical MC strategy, in green 
the antithetic variable strategy; reproduced from [9]). 


A more precise estimate of the efficiency of the approach is given on Figure 3, in which we 
have plotted the variance ratio with respect to the size of the computational domain. We see 
that the gain is not very sensitive to this size, and is at least of about 6 on this example. This 
means that, given a computational cost, the approach improves the accuracy by a factor y/6 « 
2.45. Equivalently, for a given accuracy, the computational cost is reduced by a factor 6. 



Figure 3. Efficiency of the variance reduction (same CPU time, variance ratio). 
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Our numerical results (see [6,9] for comprehensive details) show that the technique may be 
applied to a large variety of situations and has proved efficient whatever the output considered. 
In addition, we have shown in [18] that this technique carries over to nonlinear stochastic 
homogenization problems, when the problem at hand is formulated as a variational convex 
problem. In all the test cases we have considered, variance is systematically reduced. We observed 
however that the ratio of reduction is not spectacular. This has motivated the consideration 
of alternative techniques, expected to be (and indeed observed to be) more efficient than the 
antithetic variables technique. 



3. Control variate technique 

The control variate approach is a variance reduction technique known to be potentially much 
more efficient than the antithetic variable technique. It however asks to have beforehand a 
better information on the random quantity of interesf fhat is simulated. In the context of 
homogenization, the works [17,19] present a first possible investigation of fhe efficiency of fhis 
technique. 

The specific setting considered as control variate is a periodic setting slightly perturbed using 
a random field modeled by a Bernoulli variable which we now briefly describe, before turning to 
the variance reduction technique itself. 


(a) Our specific choice of control variate: a perturbation approach 

One approach, described in full details in [2-4], addressing the random material as a small 
perturbation of a periodic material, consists in considering 

Arjix^ id) = A.per {x) -t- hrj{x^ U})Cper {x), (3.1) 

where, with evident notation, Aper is a Z”*-periodic matrix modeling the unperturbed material 
and Cper is a -periodic matrix modeling the perturbation. We take 

bp{x,uj)= lQ+k{x)B^{uj), 

where the Bp are, say, independent identically distributed scalar random variables. One 
particularly interesting case (see [ 2 ^] for other cases) is when the common law of the Bp is 
assumed to be a Bernoulli law of (presumably small) parameter r): 

P(B^' = 1)=^, P(B^' = 0) = 1-77. 

A formal approach introduced in the above works (which has subsequently been studied and 
proved correct in [ 10 , 20 ]) to efficiently perform homogenization in that context starts with 
observing that, in the corrector problem 

- div [Ap {y, cj) (p -I- Vwp{y, a;))] = 0, (3.2) 

the only source of randomness comes from the coefficient Ap {y,u)). Therefore, if one knows the 
law of fhis coefficient, one knows the law of the corrector function Wp{y,u)) and therefore may 
compute the homogenized coefficient Ap, the latter being a function of this law. When the law 
of Ap has an expansion in terms of a small coefficient, so has the law of Wp. Consequently, 
Ap can be obtained as an expansion. Heuristically, on the cube Qm = [ 0 , and at order 1 
in 77 , the probability to get the perfect periodic material (entirely modeled by the matrix Aper) 
is (1 — rf)^ « 1 — N‘^r] + 0{rf'), while the probability to obtain the unperturbed material on 
all cells except one (where Ap = Aper + Cper) is — + 0{r]^). All other 

configurations, with two or more cells perturbed, yield contributions of order higher than or equal 
to rj^. This gives the intuition (and this intuition can be turned into a mathematical proof) that the 
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first order correction indeed comes from the difference between the material perfectly periodic 
except on one cell and the perfect material itself: 


+ o(r?), 

where Ap^r is the homogenized matrix for the unperturbed periodic material and 


(3.3) 


AT—>+oo 


with 


A 


1,*,N = 


j^(^per + l(3Cper)(Vu)j + Ci) — Aperi^Wi + ei)j , (3-4) 


where is the corrector for Aper (i.e. the solution to (1.1)), and wf^ solves 

-div ^(,4per + lQC'per)(V«;f^ + ei)j = 0 in Qjq, is Qjv-periodic. 


The approach has been extensively tested. It is observed that, using the perturbative approach, 
the large N limit is already very well approached for small values of N. The computational 
efficiency of the approach is clear: solving the two periodic problems with coefficients Aper and 
,4per + IgCper for a limited size N is much less expensive than solving the original, random 
corrector problem for a much larger size N. 

When the second order term is needed, configurations with two defects have to be computed. 
They all can be seen as a family of PDEs, parameterized by the geometrical location of the 
defects. Reduced basis techniques have been shown in [15] to allow for a definite speed-up in 
the computation. 


(b) Variance reduction 


We now again consider the setting defined by (3.1), except that, now, the parameter p of the 
Bernoulli law is not taken small. The expansion technique employed in Section (a) is therefore 
inaccurate. It can however serve for the construction of a control variate, useful to reduce the 
variance. 

Determining the field A{x,uj), given by (3.1), on the trimcated domain Qj^ amounts to 
drawing Bp{oj) in each cell Q + k in Qjq. This allows to compute the associated (approximate) 
homogenized coefficient A*jq (uj) from the solution to the corrector problem (3.2) truncated on Qjy. 
In parallel to this task, we reconstruct from the specific realization of the set of (ca) a field that 
is used as a control variate. More precisely, we set 

C*n{uj) = A*n{uj) - p [A;er + - E \A*per + At’^(w)]) . (3.5) 

In this formula, 

a*,N , S_ 1 Tjk, N ,ldef 

k+Q(1Qn 

where is fhe determinisfic coefficienf corresponding to the case of one defect located 

at position fc in Qjv (it is actually independent of k and equal fo 2li * defined by (3.4)). 
The parameter p in (3.5) is a deterministic parameter, a classical ingredient of control variate 
techniques, which is optimized in terms of fhe esfimated variances of fhe objecfs at play. It is 
crucial to note that the expectation of A*'^(uj) is analytically computable. Since by construction 
E (C)^) = E {A*j^), the technique then consists in approximating the former (thus the latter) by an 
empirical mean. The theoretical study and the numerical tests in [17] show that the variance of 
is smaller than that of A'^, and hence fhaf the quality of the approximation is improved. 

As an illustration, we use a similar case as in Section 2, namely (2.1) with a = 3 and /3 = 23. 
This case falls within the framework (3.1) with rj = 1/2. This is hence not a perturbative setting. 
Applying the above strategy based on (3.5) provides the results of Figure 4, where the variance is 
reduced by a factor close to 6, that is, comparable to the technique of anfifhetic variables. 
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Figure 4. Estimation of A*^ together with its confidence intervai (computed using M = 100 i.i.d. realizations), for the 
ciassicai MC simuiation (in biue) and with the controi variate approach (3.5) (in black) (reproduced from [17]). 


It is also possible to use a second order expansion with respect to rj in (3.3), and include in 
the control variate both terms, namely the deterministic coefficients corresponding to the case of 
one and two defects in Qjv- Here, one needs additional parameters playing the role of p above, in 
order to ensure substantial variance reduction (see the details in [17]). The variance reduction of 
such a case, of the order of 40, is represented on Figure 5. 



Figure 5. Estimation of together with its confidence intervai (computed using M = 100 i.i.d. realizations), for the 
classical MC simulation (in blue) and with the second-order control variate approach (in red) (reproduced from [17]). 
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4. Special Quasi-Random Structures 


The variance reduction approach we now overview has been originally introduced by other 
authors for a slightly different purpose in atomistic solid-state science [21-23]. It carries the 
name SQS, abbreviation of Special Quasirandom Structures. The approach has been adapted to the 
homogenization context in [16,19] to which we refer the reader for a more detailed presentation. 

(a) Motivation and formal derivation of SQS conditions 

In order to convey to the reader the intuition of the original approach, we first consider here 
a simple one-dimensional setting, which illustrates the difficulties of a generic problem. We 
consider a linear chain of atomistic sites of two species A and B that interact by a nearest- 
neighbour interaction potential Vaa, Vab and Vbb- 

In order to compute the energy per unit particle of that atomistic system, one has to consider 
all possible such infinite sequences, and for each of them its normalized energy 



(4.1) 


where Xi denotes the species present at the i-th site for that particular configuration. The 
"energy" of the system is then defined as the expectation of (4.1) over all possible configurations. 
The approach introduced in [21-23] consists in selecting specific truncated configurations 
(Xj)_jv<i<Ar of atomic sites that satisfy statistical properties usually obtained only in the limit of 
infinitely large N. 

The first such statistical property is the volume fraction, namely the proportion of 
species (A,B) present on average: one only considers truncated sequences (2fi)-Af<i<Af that 
exactly reproduce that volume fraction. Similarly, one may only consider tnmcated sequences 
(Xi)_ 7 v<i< 7 v that, in addition to exhibiting the exact volume fraction, have an average energy 



1 


2N + 1 , 


i=-N 


quantities of interesf. 

Mathematically, this selection of suitable configurations among all possible configurations 
amounts to replacing the computation of an expectation by that of a conditional expectation. 

The above simplistic model can of course be replaced by more elaborate models, with more 
sophisticated quantities to compute, and more demanding statistical quantities to condition 
the computations with. The bottom line of the approach remains the same, and we now 
describe its adaptation so as to construct a variance reduction approach for numerical random 
homogenization. 

To start with, we assume that the matrix valued random coefficient A present in (0.1) reads as 


Arj{x,uj) = Co(x,oj) + r/ x(x,cj)Ci(x,cj) 


(4.2) 


for some presumably small scalar coefficient r/, and where we assume that Cq and Ci are 
two stationary, coercive, uniformly bounded matrix fields, that Cq — Ci is coercive, and that 
X is a stationary scalar field with values in [—1,1]. Under these assumptions the matrix Arj is 
stationary, boimded and coercive, uniformly with respect to uj. Since rj is small, Arj is intuitively 
a perturbation of the matrix valued field Cq{x, uj). 

As already performed above, we may expand all quantities of homogenization theory 
in powers of the small parameter p. In particular, the approximations Vw^ and A^’^ of, 
respectively, the corrector Vwr/ and the homogenized matrix A^ on the truncated domain Q^, 
can be expanded in powers of ip. 


Vw^{-,uj) = VWo {■,uj) + rjVui {-, 1 . 0 ) + ri'^Vu 2 {■,uj) + o{rj‘^) 

= Al’^{uj) + vAl’^{uj) + r,^Al’^{uj) + oir,^). 


(4.3) 
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Inserting these two expansions in (1.5) and (1.4), one easily sees that 

— div C'o(p + Vwj^) = 0 inQjVi 'Wq is Qjv'Periodic, 

— divCoVu^j^ = div |^xCi(p + is Q^r-periodic, 

— divCoVu^ = div inQjVi is Q^r-periodic, 

and that the random variables Aq’^( uj), (oj) and A^’^(oj) read as 


Aq’ {u)p = 


A-k^N / \ 

{u)p = 


A-k^N / \ 

^2 (^)P = 


IQatI 

1 

IQatI 

1 

IQatI 


Co{-,uj){p + Vwo 


xi-,uj)Ci{-,uj){p + \7wo ^ 


IQatI 


Co(-,a;)Vuf(-,^), 


Q^ 


xi-,uj)Ci{-,oj)\/ui {■,UJ) + ^ 


IQtvI 


Co{;Oj)Vu^i;Oj). 


In line with the motivation we have mentioned above in the context of solid state science, we 
are now in position to introduce the conditions that we use to select particular configurations of 
the environment within Qj^. 

For finite fixed N, we say that a configuration lu satisfies the SQS conditions of order up to k if, 
for any 0< j <k, the coefficient A*'^ {uj) of the expansion (4.3) exactly matches the corresponding 
coefficient A* of the analogous expansion of the exact homogenized matrix coefficient A!^. More 
explicitly, we speak about the SQS condition of 


order 0 if A^'^(lo) = Aq, that is to say, for any p G R'*, 


1 


IQatI 


order 1 if A^’^(a;) = A^, that is to say, for any p G R'*, 
1 


Co{x,uj){p + Vwq {x,uj))dx = E 


Q? 


Uq 


Co(p +Vwo) 


(4.4) 


IQtvI 


Qi' 


{x{x,^o)Ci{x,uj){p + Vwq {x,uj)) + Cq{x,uj)Vui {x,w)^ dx 

xC'i(p + Vwo) + CoVui 


= E 


Q 


(4.5) 


order 2 if A^'^(cu) = ^4^, that is to say, for any p G R'*, 
1 


IQtvI 


Qi' 


(^xi^!'^)Ciix,uj)\/u^(x,uj) + Co{x,lu)\/u 2 (x,uj)^ dx 

= E xCiVui + C'oVu2 


(4.6) 


It is easily observed that using such particular configurations that satisfy the SQS conditions 
of order up to k we have, in the perturbative setting considered here, 

A*r,’^{u:)-A^^ = o{ri>^). (4.7) 

Taking the expectation over such configurations therefore formally provides a more accurate 
approximation of A^. Qf course, the purpose is to apply the approach beyond the perturbative 
setting. A property such as (4.7) cannot be expected any longer since the homogenized matrix 
A* is no longer a series in a small coefficient that encodes a perturbation. Nevertheless, it 
can be expected that selecting the configurations using these conditions may improve the 
approximation, in particular by reducing the variance. 
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To make the computation of the right-hand sides of the above conditions practical (since in 
theory they can only be determined using an asymptotic limit, and are therefore as challenging 
to compute in practice as A* itself), we restricf the generality of our setting. We assume that, 
in (4.2), Co (a;, a;) = Cq is a deterministic, constant matrix, Ci(x,a;) = Ci{x) is a deterministic, Z'^- 
periodic matrix, and that x{x,oj) = Xj, (cj)lfc_|_Q (x), where Xj^(oj) are identically distributed. 


not necessarily independent, bounded random variables. For the sake of simplicity, we also 
assume here that 

E[Xol=0 

and refer to [16,19] for more general cases. After a tedious but not complicated calculation (the 
detail of which is provided in [16,19]), we obtain that the two conditions (4.5)-(4.6) rewrite as 


IQwl 


Xfe(a;) — 0, 


feGZ'^nQw 


\Qn\ 


^ nXoXk]!^ 


(4.8) 

(4.9) 


fe.feQwnZ'i 


fceZ'i 


respectively, where 


CiV^iand/f, = 


fc+Q 


Ci{x)X(f>^ {x — fc)dx. In these expressions. 


Q+j 


01 is the (imique up to the addition of a constant) solution in VxG (l2(E‘^))‘^| 

to 


- div [CoV0i] = div [igCip] ir 


while (j>^ is the (unique up to the addition of a constant) solution to 

- div ^CoX4>i ] = div [igCip] inQiVi 0]^ is Qjv-periodic. 


The conditions (4.8)-(4.9) are called the SQS 1 and SQS 2 conditions. On the other hand, in the 
parhcular setting chosen, condition (4.4) (SQS 0, in some sense) is easily seen to be systematically 
satisfied when N is an integer and the truncated approximation of (1.3) that is chosen is the 
periodic approximation (1.5). 


(b) Selection Monte Carlo sampling 


The classical Monte Carlo sampling consists in successively generating a random configuration ujm, 
solving the truncated corrector problem (1.5) for that configuration, computing A*fq{uim), and 

1 ^ 

finally computing the empirical mean := — A*j^{u}m) as an approximation for A*. 


m—1 

In our selection Monte Carlo sampling, we systematically test whether the generated 
configuration satisfies the required SQS conditions, up to a certain tolerance, and reject it if it 
does not, before solving the corrector problem (1.5) for that configuration and letting it contribute 
to the empirical mean. 

In full generality, the cost of Monte Carlo approaches is usually dominated by the cost of 
draws, and therefore selection algorithms are targeted to reject as few draws as possible. In 
contrast, in the present context where boundary value problems such as (1.5) are to be solved 
repeatedly, the cost of draws for the configuration is negligible compared to the cost of the 
solution procedure for such boimdary value problems. Likewise, evaluating the quantities present 
in (4.8)-(4.9) is inexpensive. Therefore, fhe purpose of the selection mechanism is to limit the 
number of boundary value problems to be solved, even though this comes at the (tiny) price of 
rejecting many configurations. We also note that, as for any selection procedure, our selection may 
introduce a bias (i.e. a modification of the systematic error in (1.6)). The point is to ensure that the 
gain in variance dominates the bias introduced by the selection approach. 

We have studied the approach theoretically in [16,19]. It is shown therein that the estimator 
provided (at least the simplest variant of our approach) converges towards the homogenized 
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coefficient A* when the truncated domain converges to the whole space. The efficiency of the 
approach is also theoretically demonstrated for some particular and simple situations (such as 
the one-dimensional setting). A comprehensive experimental study of the approach has been 
completed. In particular, since it is often necessary to enforce the desired conditions only up 
to some tolerance, we have investigated in [16,19] how this tolerance affects the quality of the 
approximation and the efficiency of the approach. We have observed that the approach is robust 
in this respect. 

We include here a typical illustration of the efficiency of the approach. We again use a similar 
case as in Section 2, namely (2.1) with a = 1/2 and /3 = 3/2. Considering only configurations 
that exactly satisfy (4.8), we obtain the results shown on Figure 6. It is also possible, among the 
configurations that exactly satisfy (4.8), to select configurations that satisfy as best as possible the 
condition (4.9). In practice, we generate 2000 configurations that exactly satisfy (4.8) and select 
among them the 100 configurations for which the difference between the left and the right-hand 
sides of (4.9) is the smallest. We then obtain the results shown on Figure 7. 



50 


Figure 6. Estimation of AJj together with its confidence interval (computed using M = 100 i.i.d. realizations) as a 
function of N, for the classical MG simulation (in black) and with the SQS approach based on (4.8) (in red) (reproduced 
from [16]). 


In the case considered here (for which the contrast in the field A is equal to 3), the variance 
is reduced by a factor 20 when using configurations that exactly satisfy (4.8), and by a factor 300 
if (4.9) is enforced as well. To compare this variance reduction approach with the two previous 
ones, it is however needed to consider a case for which the contrast in A is similar. In that case, 
the variance is reduced by a factor of 9 when using configurations that exactly satisfy (4.8), and 
by a factor of 60 if (4.9) is enforced as well. 

In all the test cases we have considered (see [16,19] for details), we have observed that the 
systematic error is kept approximately constant by the approach (it might even be reduced), while 
the variance is reduced by several orders of magnitude. Such an efficiency is achieved at almost 
no additional cost with respect to the classical Monte Carlo algorithm. 


5. Related issues and Further research 


The studies we have reviewed above on different variance reduction approaches definitely show 
that such approaches may be very beneficial in the context of random homogenization, improving 
the accuracy while essentially preserving the computational cost. Their efficiency, measured as the 
actual ratio between the variance of a quantity computed with a direct Monte-Carlo approach and 
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Figure 7. Estimation of together with its confidence interval (computed using M = 100 i.i.d. realizations) as a 
function of N, for the classical MC simulation (in black) and with the SQS approach based on (4.8) and (4.9) (in blue) 
(reproduced from [16]). 


that of the same quantity computed using the variance reduction approaches, varies, depending 
upon the amount of information that one has on the problem and that one inserts into the specific 
variance reduction approach. The antithetic variable approach, a quite generic approach that 
can be put in action almost without any prior knowledge on the problem considered, already 
reduces the variance by one order of magnitude, say, in the best case scenarios. Control variate 
and Special QuasiRandom Structures, both approaches that require exploiting some information 
on the problem, perform much better. Their efficiency may typically be one order of magnitude 
larger. 

Of course, the efficiency of all approaches is sensitive to the contrast present in the original 
multiscale problem. In a schematic manner, one may say that the efficiency is, approximately, 
inversely proportional to the contrast. It is an issue, since practically relevant multiscale problems 
may present a high contrast. Fortunately, there is room for improvement in the approaches and 
several ideas, some of them already explored in other contexts of the engineering sciences, some 
of them not, have not been pursued yet. 

Among possible tracks for further research, we wish to cite a couple of alternate control variate 
approaches. 

A first possible track consists in considering nonlinear convex stochastic homogenization 
problems (as those considered in [18]), and use a corresponding linear problem either as a control 
variate (in the spirit of the approaches presented in Section 3) or as a way to select particular 
configurations (as in Section 4). We do not detail here the precise construction of this linear 
model, but rather focus on how to use it in practice. Let ^ G i->- VL*(^) G R be the homogenized 
energy density of the nonlinear stochastic homogenization problem, and ^ be its 

approximation computed by considering the nonlinear cell problem on the bounded domain Qjv- 
Let A^(u!) be the homogenized matrix of the corresponding linear problem. Our aim is to use 
as a control variate for ^). Note however that we do not know the expectation 

of and hence we cannot directly use a Monte Carlo algorithm on the random variable 

W^f{uj,0 - p(c^A]V(^)^ - E ). 

However, computing is expected to be less expensive than computing W^{uj, ^), because 

the corrector problem in the former case is linear, whereas it is nonlinear in the latter case. A 
natural idea is thus to replace, in the above relation, E by an empirical mean. This 
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leads to approximate E 5)] by a mean of the form 

M M 

m—1 m—1 

where M, A4 (that we expect to be much larger than M) and p are chosen to minimize the variance 
of the approximation for a given computational cost. 

A second track for further research is to use the so-called bounds, that are routinely employed 
in Mechanics, in order to build a control variate approach. Given the computational cost for 
obtaining approximations of A*, practitioners indeed sometimes choose to avoid computing the 
actual homogenized coefficients (by solving (1.4)-(1.5)) and concentrate on bounds (namely the 
Reuss, Voigt, Hashin-Shtrikman bounds,...) on the homogenized matrix A*. 

For the sake of illustration, let us briefly review the derivation of the so-called Voigt bound. We 
assume that the random coefficient A is a symmetric matrix. This assumption is critically used in 
what follows, and more generally in the derivation of many bounds. Under this assumption, the 
matrix defined by (1.4), satisfies, for any p, 


(cj)p = inf 


{p + 'Vv)'^A{-,uj){p + Vv), veHpeAQN) 


. \Qn\ Jq„ 

and hence, by choosing r; = 0 in the above problem, we obtain that 

1 


^Ar(w) < 


IQwl 


Qn 


A{-,uj). 


The average of A(-, u) over hence provides an upper bound on A'^(uj), which is the so-called 
Voigt bound. 

In the specific case of two-phase composite materials (made of two phases denoted A and B), 
where the random coefficient is given, with obvious notations, by 


A{x, oj) = x{x,uj) A + {1 - x{x,oj)) B, 


where x is the characteristic fimction of the phase A, more elaborate boimds have been proposed, 
including the so-called Hashin-Shtrikman boimds. We refer e.g. to [1] for more details. The idea 
we are currently pursuing is to use these bounds not as an approximation for A^(a;), but as a 
control variate. 
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